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Linearly forced isotropic turbulence 

By T. S. Lundgren f 


Stationary isotropic turbulence is often studied numerically by adding a forcing term 
to the Navier-Stokes equation. This is usually done for the purpose of achieving higher 
Reynolds number and longer statistics than is possible for isotropic decaying turbulence. 
It is generally accepted that forcing the Navier-Stokes equation at low wave number does 
not influence the small scale statistics of the flow provided that there is wide separation 
between the largest and smallest scales. It will be shown, however, that the spectral 
width of the forcing has a noticeable effect on inertial range statistics. A case will be 
made here for using a broader form of forcing in order to compare computed isotropic 
stationary turbulence with (decaying) grid turbulence. It is shown that using a forcing 
function which is directly proportional to the velocity has physical meaning and gives 
results which are closer to both homogeneous and non-homogeneous turbulence. 

Section 1 presents a four part series of motivations for linear forcing. Section 2 puts 
linear forcing to a numerical test with a pseudospectral computation. 


1. Motivation for Linear Forcing 

1.1. Linearity of energy production in non-homogeneous turbulence 
In shearflow turbulence the equation for the fluctuating part of the velocity, u', is 


(9u/ 

— - + u • Vu' + u' • Vu + iT • Vu' - V • u'u' 
at 


—'S7p'/p + j'V 2 u'. 


(1) 


The third term on the left, u' • Vu appears in the turbulent energy equation as the energy 
production term, < u' • Vu • u' >. (Both angle brackets and overbars are used to denote 
averages.) In (1) it appears as a forcing term proportional to u'. This suggests that for 
isotropic homogeneous turbulence it might be appropriate to force a stationary flow with 
a driving term proportional to the velocity. Of course, for isotropic turbulence there is 
no mean velocity gradient, so the only way to have the flow be isotropic and stationary 
is to have the forcing be isotropic. It is proposed here to use 

— + u • Vu = -S/p/p + iA7 2 u+ f (2) 


with the driving force 


f = Qu , 


( 3 ) 


where Q is a constant. The prime on the velocity has been omitted here and henceforth 
with the understanding that < u >= 0. The turbulent energy equation is now 


1 d < u • u > 

2 dt 


—e + Q < u • u >, 


( 4 ) 


where 


e = —v < u • V 2 u > 


I University of Minnesota 


( 5 ) 
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is the mean dissipation rate and the last term could be called isotropic turbulent produc- 
tion. For stationary turbulence therefore (setting the time derivative to zero) 

Q = e/3t/ 2 (6) 

where U 2 =< u • u > /3 is the mean square of one component of the velocity. The 
proposal is to numerically solve 

— + u • Vu = -Vp/p + i/V 2 u + (e/3[/ 2 )u (7) 

dt 

with the objective of comparing the statistics with those of grid flow turbulence. Equa- 
tion (7) has the property that the rest state is unstable to long waves. Therefore solutions 
cannot decay to zero but must transfer energy to shorter waves in order to dissipate en- 
ergy. Some reasons are given below for expecting computational results to be comparable 
with experiments for inertial scales of turbulence. 


1.2. Linear forcing of the Karman-Howarth equation 
A derivation of the Karman-Howarth equation, following the steps given in detail by 
Landau and Lifshitz (1959), but including a forcing term as in (2), gives 

dU 2 1 dB 2 11 9 4 
dt 2 dt 6 r 4 dr r 3 ~ 


1 f) 1 C r 

- + -g / r 2 < u(x 2 , t) • f(xi, f) + u(xi, t) ■ f(x 2 , t) > dr (8) 
m or or m 

where B 2 {r, t) and B%{r,t) are second and third order longitudinal structure functions, 
X 2 = xi + r and r = |r|. With f given by (3), the last term can be written 

2 Qr 2 Ru(r,t)dr (9) 

where Rij(r,t) =< Uifx.\,t)ujfx.i + r, t) > is the velocity correlation tensor. Since the 
trace of the correlation tensor is related to the second order structure function by 

Ru = 3U 2 -^^-r 3 B 2 (10) 

2 r z or 

(9) can be integrated to the form 

2 QU 2 ~ QB 2 (r,t) . (11) 



Using this with Q given by (6), and dropping the time derivatives for stationary turbu- 
lence, gives the result 



1 <9 1 d dB 2 

6 r 4 dr 3 r 4 dr dr 


(12) 


in which e is a constant given by (5). This differs from the standard decaying version of 
the Karman-Howarth equation which is 


2 1 dB 2 11 9 4n Id 4 dB 2 

3 2 dt 6 r 4 <9r r 4 dr dr 


(13) 


where here e is a function of time given by 

. , 3 dU 2 

'« = -2 nr 


15v d 2 B 2 

2 dr 2 


o 


(14) 
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The Kolmogorov r 2 / 3 law may be derived from (12), by the method of matched asymp- 
totic expansions, in a manner similar to that employed by Lundgren (2002) to get that 
law from (13) but without the necessity of using similarity in time. The result is the 
same: 

B 2 = C 2 U 2 (r/L) 2 ' 3 = C 2 (er) 2 ' 3 (15) 

where L = U 3 /e and e and U 2 are independent of time now. 

Equation(12) can be integrated to obtain B 3 in terms of B 2 , as 


B 3 = 6u 


dB 2 

dr 



2e 1 
IP P 



Using (15) for B 2 gives the simple result 


4C 2 / r \ 

-i/ 3 4 / r 

a ec 2 /, 

r \ 5/3 

R l \LJ 

“5 \l 

J + 17 L 

l) 


(16) 


(17) 


This may be rewritten using the Taylor microscale A to scale r: 


B 3 = 



(18) 


Here Rl = UL/v\ Rl = R\/ 15 relates R\ and Rl and \/L = 15 /R\ relates A and 
L. The corresponding result for decaying turbulence (Lundgren 2002; Linclborg 1999) is 
almost the same: 


B 3 = 



1 -R 


-2/3 

A 




(19) 


where n is the energy decay exponent ( U 2 a t~ n ). Equations (18) and (19) would be 
exactly the same if n = 2, which is not a realistic decay exponent; n = 1.2 is often 
observed and n = 4/3 is the maximum value possible. These equations give a Reynolds 
number correction to the Kolmogorov “4/5” law showing that it is approached slowly as 
R\ — » oo. When C 2 = 2 and n = 1.2 the compensated forms have maxima at r / A = 1.23 
(for linearly forced turbulence) and r/A = 1.11 (for decaying turbulence). 

The similarity of these equations can be seen in Figure 1 where they are plotted for 
R\ = 350 and compared with an experiment in a turbulent jet. 


1.3. Normalized decaying turbulence 

Consider (7) again, with and without the isotropic forcing term. In the stationary case 
with forcing introduce dimensionless variables 

v = u /U, x = r/L, r = tU/L (20) 

where U and L = U 3 /e are constants. A simple change of variables gives 

-^+v-Vv = -VP + i?7 1 V 2 u+-v . (21) 

dr L 3 

Now consider the case without the forcing term, with the change of variables 

v = u/U(t), x = r /L(t), t= f U/Ldt (22) 

J o 
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Figure 1. Compensated third order structure function vs. r/ A at R \ = 350. The data points are 
from a laboratory jet (Gagne, 2002) at the same value of R\. The long dashed curve is from (19). 
The dashed curve is the linear forcing result from (18). The solid line is from the computation 
of section 2 using the R scaling law to extrapolate from R\ = 170 to R\ = 350. 

The equation transforms to 

^ + vVv = -VP + i^ 1 V 2 v-|§v+^x-Vv. (23) 

or U z U 

Using e = — § and L = U 3 /e the coefficient of v is — so 

o 1 L 

-^ + V • Vv = -VP + i?7 1 V 2 v + -V +— x-Vv. (24) 

or 3 U 

While this is not exactly the same as (21) because of the last term, it is apparent that 
energy decay has the effect of an isotropic forcing term. Note also that if U 2 oc t~ n then 
e oc t~ n ~ 1 and L = U 3 /e oc < 1-n / 2 . So L would be zero if n = 2. In this case the equations 
would be the same except for the time dependence of Rl- 

1.4. Comparison with low wavenumber forcing 
In this subsection a more general forcing function is applied to the Karman-Howarth 
equation in order to show that the forcing range has a significant effect on the third or- 
der compensated structure function and presumably on any inertial range statistics. This 
analysis is carried out in detail in Appendix I and only summarized here. A Gaussian 
filter is applied to the velocity field, filtering out a variable part of the high wavenumber 
content. This filtered velocity is used as a forcing function for the Karman-Howarth equa- 
tion. In one limit there is uniform linear forcing. The results of a calculation are shown 
in Figure 2a for different values of KL. Here I\ is the width of the filter; wavenumbers 
greater than K are filtered out. When K is small only low wavenumbers are forced, while 
if K — > oo forcing is uniform over all wavenumbers. The figure shows the compensated 
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Figure 2a. Third order compensated structure functions for variations of the forcing range com- 
puted from (A17). The curves are (from the bottom) for KL = 1000, 100, 25, 10, 5, 3. The lowest 
curve (KL = 1000) is close to the linear forcing result from (18). All the curves are computed 
for R x = 1000. 



r/rj 

Figure 2b. Solid curves, third order compensated structure function vs. r/g for KL = 5 and 
R\ = 125,284,381,460 from the bottom up. This compares with the low wavenumber forced 
computation of Gotoh et. al. (2002) at the same Reynolds numbers (compare their Fig. 12). The 
maximum values shift to the right and increase with increasing R\ in a similar manner. The 
maxima are .70, .76, .77, .78 comparing with .66, .77, .78, .76 from Gotoh. The dashed curves 
are for linear forcing at the same values of R\ from (18) showing the differences. 
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Time 


Figure 3. Time History of R\ with linear forcing. Numerical data for the following figures was 
taken over a time range of about .5 near T=10, where R\ = 170. 


third order structure function for different values of KL for R\ = 1000. When KL = 1000 
the curve is nearly the same as computed from (18). As KL decreases the maximum value 
moves to the right and increases towards 4/5, greatly distorting the structure function 
in the inertial range. 

Figure 2b shows compensated structure functions for several values of R\ for the single 
value KL = 5. This was constructed in order to compare with the high resolution DNS 
of Gotoh et. al. (2002), which was driven by white noise in a low wavenumber band 
which corresponds roughly to a Gaussian filter with width KL = 5. (In their numerics 
at R\ = 460 L ~ 2.5, making K = 2 which is in approximate accord with their forcing 
range of 1 < k < V&-) 


2. Pseudospectral Computation with Linear Forcing 

Computations have been carried out with a Rogallo/Wray pseudospectral code, mod- 
ified slightly to accomodate linear forcing. This was done with (256) 3 resolution with 
v = .003, box size 27 t on a side, and time scale set by taking Q = 1. Figure 3 shows 
the time history of R\ during a lengthy run. As can be seen the computation is not 
exactly stationary (but should be statistically stationary after an initial transient). The 
instantaneous value of e/3 U 2 (not shown here) fluctuates about unity by about ±20%. 
There are fairly stationary periods, however. Data for the following figures was taken in 
a relatively stationary period near T = 10 where R\ = 170. 

Figure 4 shows the energy spectrum versus kg. The spectrum is a little shallower than 
fc -5 / 3 (the upper straight line on the figure). Mydlarski and Warhaft (1996) indicate 
approximately fc -145 at this Reynolds number for grid flow turbulence. 
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kr) 


Figure 4. Energy spectrum in Kolmogorov variables for R\ = 170. The upper straight line is 
k ~ 5 ^ 3 and the lower one is k ~ 1A 5 . 



Figure 5. Compensated even order structure functions versus r/X for R\ = 170. Solid lines are 
curve fit using (25). The solid straight lines are an extrapolation to R\ = oo using V 2 , 14, Ve- 


in Figure 5 compensated structure functions of second, fourth, and sixth order are 
plotted versus r/X. The minimum separation between two points (27 t/ 256) is very nearly 
.2A. The output was taken by averaging over all pairs of points along the x direction with 
separation r (which is a multiple of .2A). For each structure function there are six curves 
(the lighter long-dashed curves on the figure), one each for times separated by 1000 time 
steps. The total elapsed time for 5000 time steps is about .5 time units on Figure 3. The 
spread of these curves indicates the random nature of the output even after the spatial 
averaging. 
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Figure 6. Odd order compensated structure functions for R x = 170. Curve fit same as in Figure 
5. 


The pth order compensated structure function should take the form (Lundgren 2003) 

B P /(er) p ' 3 = V p (l - R x 2/3 {A p {r/X) 2 / 3 + B p (r/ X)~ 4/3 )^j . (25) 

This gives Reynolds number corrections to Kolmogorov (1941) theory and can be re- 
garded as a scaling law. If the parameters V pi A p , B p were determined for some Reynolds 
number (numerically or by experiment) one could extrapolate the data to another Reynolds 
number. The coefficients for p = 2, 3, 4, 5, 6 were determined as follows. The six time sets 
were truncated to about .5 < r/A < 5 in order to get the upper parts of the sets, roughly 
centered about the maxima. These were accumulated into into a single set, which looks 
like a lot of scattered experimental points. Then using a nonlinear curve fitting algorithm 
(on xmgr) the coefficients V p ,A p ,B p were determined for R x = 170 ( R x 2 '' 3 = .033). The 
values are shown in the table below. 


p 

V P 

A p 

B p 

(r / A) max 

V p /pl 

2 

2.11 

2.62 

6.15 

2.17 

1.05 

3 

-.824 

6.56 

3.75 

1.07 

-.137 

4 

15.24 

4.85 

7.65 

1.78 

.635 

5 

-16.91 

7.81 

4.53 

1.08 

-.141 

6 

197.5 

6.56 

6.87 

1.45 

.274 


Table 1. Coefficients for (25). Fifth column is position of maximum of B^/er. Sixth 
column compares V p with p\. 


The fit curves are shown as heavy lines which approximate the data fairly well. These 
should be thought of as time averages of the six curves of each set. The horizontal straight 
lines above each set are the curves B p /(er ) p / 3 = V p , representing an extrapolation to 
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Figure 7. Second and third order structure functions repeated from figures 5 and 6 on a semilog 
plot in order to better show the nature of the curve fit. 

R\ = oo. In particular note that the horizontal line above the second order stucture 
function, V 2 = 2.11, is an acceptable value for that asymptote. 

The third and fifth order structure functions were treated in the same way and shown 
in Figure 6. The individual sets show much more scatter than the even order structure 
functions because of cancellation between the left and right sides of the almost symmetric 
pdfs. Note the value V 3 = .824. This should be exactly .8, the Kolmogorov “4/5” law. 
Figure 7 repeats the B 2 and B 3 curves from Figures 5 and 6 on a semilog plot. 

The maxima of the compensated structure functions have positions which scale with A. 
For even orders the positions shift towards smaller values with increasing order. For odd 
orders there is not much shift. The rapid increase of the magnitude of the even orders 
with increasing order required that they be presented on a log-log plot in order to show 
them on the same figure. A rapid increase of moments like p\, as seen in the table, is a 
signature for exponential tails on the pdf. This is observed approximately in Figure 8. 

An example of extrapolation for the third order structure function is shown in Figure 
1, where the curve fit values of V 3 , A 3 , B 3 are used to extrapolate to R\ = 350 ( R = 
. 020 ). 

Figure 8 is the final figure. This shows the velocity pdf versus v/a where v is the velocity 
difference between two points along the x direction and a is the standard deviation, 
a =< v 2 > 1 / 2 . In these variables the average and the second moment are both unity. 
The outermost curve is as close to the pdf of the velocity derivative as can be reached 
with this resolution. The fluctuations in the tails were ameliorated to some extent by 
accumulating the pdfs over 10 successive time steps (in addition to spatial averaging). 


3. Conclusions 

There were two objectives to this research. The first was to show that it is desirable and 
useful to compute stationary isotropic turbulence with a forcing function which is a simple 
linear function of velocity, which thus forces uniformly at all wavenumbers. This was done 
by analysis of of a forced Karman-Howarth equation, which allows calculation of third 
order structure functions in terms of second order structure functions. This allowed a 
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Figure 8. pdf of velocity difference vs. via, a =< v 2 > 1,/2 The curves are for values of the 
separation. From the outside r/\ = 2., .4, .6, 1.0, 1.4, 2.0, 2.4 


favorable comparison between decaying isotropic turbulence and linearly forced isotropic 
turbulence in subsection 1.2. 

By applying filtered linear forcing to the Karman-Howarth equation it was shown that 
forcing at low wavenumber has a large effect on third-order structure functions, and would 
very likely influence all inertial range statistics, in partcular it could affect the computa- 
tion of anomalous exponents. The third order structure function with low wavenumber 
filtered forcing was compared with high quality DNS (Gotoh et. al., 2002) and linear 
forcing was compared with experiments of Gagne (2002). The differences between low 
wave number forcing and linear forcing were considerable, as seen in Figure 2b, with low 
wavenumber forcing approaching the four-fifths law too rapidly with increasing Reynolds 
number. 

The second objective was to compute a moderate resolution DNS of box turbulence 
with linear forcing to show the feasibility of such a computation. This was done in section 
2. Compensated structure functions of second through sixth order were computed at 
R\ = 170. It was shown that these could be fit with an R x scaling law, determining 
the coefficients to fit equation (25). (The analysis which produced this law (Lundgren 
2003) applies equally to the forced Navier-Stokes equation.) The usefulness of the scaling 
law was demonstrated in Figure 1, where the third-order structure function computed at 
R\ = 170 was extrapolated to R\ = 350 (the solid curve in Figure 8 was extrapolated 
to the solid curve in Figure 1), where it compares with the Karman-Howarth results and 
with experiment. 
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5. Appendix I. Filtered Forcing. 

With a general forcing term, the integrated Karman-Howarth equation is 

B 3 = - ^er - ^ I rdr f < u 2 • fq + ii* • f 2 > r 2 dr . (Al) 

dr 5 r 4 J 0 J 0 

It is desired to develop a general case which includes forcing at small wavenumbers 
and also broader forcing. A simple model takes the forcing function proportional to the 
velocity, but with a variable part of the high wavenumber end filtered out with a Gaussian 
filter. A filtered velocity in Fourier space may be written 

u < (k) = g K (k)u(k) (A2) 

where the filter and its transform are 

9k(x) = / exp(ik • x)gx(k)dk (A3) 


(2t r) 3 

It is assumed that gx{ 0) = 1 so that 


9 R {k) = J exp(— ik • x)g K (x)dx 


J g K (x)dx = (2tt) 3 


(A4) 

(A5) 


For a Gaussian filter 

g K = exp(— ,5fc 2 /A" 2 ) . (A6) 

This filters out wavenumbers k > K for 0<AT<oo;A' = oois the uniform forcing case. 
The transform of gx is 

9 k{x) = (27t) 3 / 2 AT 3 exp(— K 2 x 2 /2) . (A7) 

In general the filtered physical space velocity is given by a convolution: 


u<(x) = 


1 

(2tt) : 


J 5 ic(| x -s|)u(s)cis . 


The forcing function will be assumed to be f = Qu < , where, since e =< u • f > 

Q = e/ < u • > . 


(A8) 

(A9) 


From (A8) 

< u • u< >= J 9k( l x - s l) < u ( s ) ’ u ( x ) > ds (^ 10 ) 

where 


< u(s) • u(x) >= Ra( \x - s|) (All) 

is the trace of the correlation function given by (10). By using (A9) and (A10) the force 
correlation in (Al) may be written 


< u 2 • fi + ui • f 2 >= 2e 


/ gx{\x - s\)Ru(s)ds 
J g K (s)Ru(s)ds 


(A12) 


B3 = 6u 


dlh 

dr 


Substituting this into (Al) gives 



(A13) 
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|er(27r) 3 3t/ 2 + / 0 ' rdr / 0 ' r 2 dr f g K (\s — r|) ( — §4, £s 3 B 2 (s))ds 

(2tt) 3 3U 2 + f g K (s)( - ±££s 3 B 2 (s))ds 

where the integration over the 3 U 2 part of Ru was done using (A5) and rdr J ( ' r 2 dr = 
r 5 / 15. In the upper integral do the integration over the angle variables with the Gaussian 
filter, obtaining 


where 


B(s,r) = 


gx{ |s — r|) sin {6)d6d(j> = (2n) 5 ^ 2 KB(s, r) 

( ( exp(— .5A' 2 (s — r) 2 ) — exp(— ,5A' 2 (s + r) 2 )) 
l sr 


(AU) 


(A15) 


is a temporary notation. Now do the integration by parts on s, which shifts the differen- 
tiation onto B(s,r). The following identity may be proved 


dB(s, r) 


1 dA(s,r) 


(- 4 . 16 ) 


ds K 2 s 2 r 2 dr 

where 

A(s,r) = (1 — K 2 sr) exp(— .5 I\ 2 {s — r) 2 ) — (1 + K 2 sr) exp(— .5/v 2 (s + r) 2 ) . 

Because of the 1/r 2 factor in (A 16) the inner f r 2 dr may be carried out, resulting in the 
near-final form 

B 3 _ 6L 2 db 2 | + (27r) -1//2 /\ -1 ^ / 0 r rdr f 0 °° A(s, r)b 2 (s)sds 
er Rl dr 


(-417) 


1 — (27 r) 1 / 2 A' 5 | J 0 °° exp(— ,5K 2 s 2 )b 2 (s)s 4 ds 

where b 2 = B 2 /U 2 . The double integrals have to be integrated numerically with an 
assumed function for B 2 . An appropriate form is 


B 2 = 2t/ 2 tanh(.5C' 2 (r/A) 2/3 ) 


(-418) 


which gives the Kolmogorov two-thirds law for small r/A and tends to the proper limit 
2 U 2 as r/L — > oo. 
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